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Abstract: The extended Heine-Stieltjes polynomials associated with a special Lipkin-Meshkov- Click (LMC) 
model corresponding to the standard two-site Bose-Hubbard model are derived based on the Stieltjes correspon- 
^ - ^ ^ dence. It is shown that there is a one-to-one correspondence between zeros of this new polynomial and solutions of 

' the Bethe ansatz equations for the LMC model. A one-dimensional classical electrostatic analogue corresponding 

to the special LMG model is established according to Stieltjes early work. It shows that any possible configuration 
of equilibrium positions of the charges in the electrostatic problem corresponds uniquely to one set of roots of the 
Bethe ansatz equations for the LMG model, and the number of possible configurations of equilibrium positions of 
the charges equals exactly to the number of energy levels in the LMC model. Some relations of sums of powers 
and inverse powers of zeros of the new polynomials related to the eigenenergies of the LMC model are derived. 
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^ ! 1 Introduction 

ff^ , It is well-knovifn from a number of studies [1-3] that the Lipkin-Meshkov-GUck (LMG) model [4-6] has a 

' rich phase structure that is related to a broad range of phenomena in a variety of physical applications, 

such as analyses of spin systems [7], Bose-Einstein condensates [8], etc. Furthermore, it has been proven 
[9,10] that the model is exactly solvable. A special case of the theory is related to the standard two- 
site Bose-Hubbard model [11] vi^ith a Hamiltonian that is equivalent to that of a paramagnet of the 
easy-axis type in a transverse magnetic field [12]. The quantum dynamics and the energy gap of the 
model were studied in [13-14]. The quantum critical behavior and a special case of the model were 
' discussed in [15-16]. As shown in [9-11, 12-16], exact solutions of the model can be determined using the 

^ . Gaudin-Richardson or Bethe ansatz method, in which eigenstates of the system are written in terms of 

a product of spectral parameter-dependent operators. The unknown spectral parameters must satisfy a 
set of coupled non-linear equations, called the Bethe ansatz equations (BAEs). Solutions of the BAEs 
simultaneously determine the eigenenergies and the corresponding eigenstates [9-11]. 

It should be noted that according to the early work of Stieltjes [17-20] there is an important one-to-one 
correspondence between every set of BAEs and a set of orthogonal polynomials. And furthermore, as 
shown by Szego [21], the Stieltjes correspondence can be used to define a large class of classical orthogonal 
polynomials. Roots of these BAEs are zeros of the corresponding polynomials, which can be interpreted 
as stable equilibrium positions for a set of free charges in an external electrostatic field. This link between 
Richardson's BCS pairing model for nuclei and the corresponding electrostatic problem was established 
in [22] based on an earlier unpublished preprint of Gaudin, which was then made clearer in [23] . A much 
more general approach to the pairing model was shown in [24] and [25]. The purpose of this paper is 
to show that the Stieltjes correspondence related to the solutions of a special LMG model corresponding 
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to the standard two-site Bose-Hubbard model gives rise to a new polynomials, which is now called the 
extended Heine-Stieltjes polynomials . 



2 Stieltjes correspondence 

There is a large class of polynomials, called Heine-Stieltjes polynomials y{z), which satisfy a second-order 
Fuchsian equation 

A{x)y"ix) + B{x)y'{x) + V{x)y{x) = 0, (1) 

where A{x) is a polynomial of degree m with A{x) = n"=i(^' ~ '^i)^ ^{^) ^ polynomial of degree m — 1 
such that for a set of real positive parameters {7i} and another set of real parameters {a,} 



B{x)/A{x)=J2-rh' (2) 

. T X Hi 
1=1 



and V{x) is an unknown, but to be determined polynomial of degree m — 2 that is allowed to depend 
on the solution y{x). The latter are often called Van Vleck polynomials [26]. The case with m = 2 
corresponds to the hypergeometric differential equation, while the case with m = 3 corresponds to the 
Heun equation [27-28]. 

Such Heine-Sieltjes polynomials and properties of their zeros have been studied extensively. If the 
polynomials A{x) and B(x) are algebraically independent, i.e., they do not satisfy any algebraic equation 
with integer coefficients, Heine proved that for every integer k there exist at most a{k) = {k + m — 
2)!/((to — 2)!fc!) different Van Vleck polynomials V{x) such that y{x) has a polynomial solution of degree 
k. As summarized in Szego's work on orthogonal polynomials [21], for every set of nonnegative integers 
(fci, ^2, • • • , krn), there are uniquely determined real values of the parameters for the Van Vleck polynomials 
V{x) such that y(x) has a polynomial solution with kj simple zeros in the open interval (aj_i, a^) for each 
j = 1, 2, • • • , m. The polynomial y{x) is uniquely determined up to a constant factor, and has the degree 
k = ki + ■ ■ ■ km- As noted by Volkmer [29], this existence and uniqueness theorem can also be derived 
from the general Klein oscillation theorem for multi-parameter eigenvalue problems shown in [30]. 

If y{x) is a polynomial of degree k with simple zeros {xi,X2, ■ ■ ■ , Xk}, one may write y{x) as 

k 

y{x) = l[{x-xj). (3) 



Thus, one can easily check that y(x) will satisfy 



y'ixi) ^ X, 



(4) 



at any zero Xi. Combining (1), (2), and (4), one obtains the following important relations among the 
zeros: 
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for i = 1, 2, • • • , A;. It should be noted that the BAEs frequently appearing in search for exact solutions 
of quantum many-body problems, such as those associated with Gaudin type systems, are similar to the 
relations shown in (5). The link between Richardson's BCS pairing model for nuclei and the corresponding 
electrostatic problem was investigated in [22] based on an earlier unpublished preprint of Gaudin, which 
was then made clearer in [23]. A much more general approach to the pairing model was shown in 
[24] and [25]. Roots of the BAEs (5) simultaneously determine the cigcncncrgies and eigenstates of the 
corresponding quantum many-body problem. As an alternative, roots of the BAEs may also be calculated 
as zeros of the corresponding polynomial y{x). In this way, a link between solutions of the Gaudin type 
for quantum many-body problems and the corresponding polynomials is established. 



3 Special LMG model and Bethe ansatz solutions 

As a simple extension of the Stieltjes correspondence, we revisit the Bethe ansatz solutions for a special 
LMG model corresponding to the standard two-site Bose-Hubbard model studied in [11, 13-15], for which 
the Hamiltonian is 

H = -t{cU + d)c) + U{c^cc'^c + d)dd)d), (6) 

where (c) and d) (d) are boson creation (annihilation) operators, the parameters t and U in the Bose- 
Hubbard model are related to the Josephson coupling and the charging energy, respectively. Following 
[11, 15], we use the unitary transformation for the boson operators with 

c=^|\{a-^h), d=^{a + ib). (7) 

Then, the Hamiltonian (6) can be rewritten in terms of a- and 6-boson operators as 

H = t{b^b - a+a) - ^US+{0)S^{0) + Uh^, (8) 

where n = a^a + b^b is the operator for the total number of bosons in the system, and 

5+(0) = 6t2 + at2, 5-(0) = 62 + a2 (9) 

are boson pairing operators. 

Let \i'i,V2) = a^'^^b^'^^\0) with j/j = or 1 for i = 1, 2 be the a- and 6-boson pairing vacuum state 
satisfying 

a'^\ui,U2) =b'^\vi,U2) =0. (10) 
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To diagonalize the Hamiltonian, we use the algebraic Bethe ansatz which implies that eigenvectors of (8) 
may be expressed as 

|n, C,'^i,'^2) = S+{x[^'>)S+{xi^^)---S+{xi^%,,,,2) (11) 
with n = 2k + 1^1 + 1/2, and 



xy-' + 1 x\^' - 1 

in which a;-^^ {i = 1, 2, • • • , fc) are spectral parameters to be determined, and C is an additional quantum 
number for distinguishing different eigenvectors with the same quantum number k. It can then be 
verified by using the corresponding eigen-equation that (11) is a solution when the spectral parameters 
x^p (z = 1, 2, • • • , fc) satisfy the following set of BAEs: 



at2 6t2 



U 



for i = 1, 2, ■ • ■ , fc, with the corresponding eigen-energy given by 

k 



E^O =2tJ2 4^' + *(^2 - i^i) + Un''. (14) 



i=l 

Hence, once the (^-fh roots {a;^^^} are obtained from Eq. (13), the eigenenergy and the corresponding 
eigenstate are thus determined according to (14), (11), and (12). 



4 The extended Heine- St ieltjes polynomials associated with so- 
lutions of the LMG model 

In order to compare to the Stieltjes correspondence, we assume t > and U ^ 0, and set 

a = ui + l/2, p = U2 + l/2, 'y = t/U. (15) 

Then, the BAEs (13) becomes 

When U > 0, all parameters a, (5 and 7 are always positive. When [/ < 0, we interchange the boson 
operator with 6^ in (12). Then, one can easily verify that such a change is equivalent to interchange 
v\ with f2 and t ^ —t m the BAEs (13), which leads to the final BAEs for the f/ < case that is the 
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same as (16) with a ^ ^ and keeping a > 0, /3 > 0, and 7 > 0. Hence, it is sufBcient to consider the 
case for ?7 > only with BAEs given by (16). 

Although the 7 = result is trivial for the LMG model, according to the Stieltjes correspondence, 
the polynomial corresponding to (16) in this case is the Jacobi polynomial P^°'~^'^~^\x) satisfying the 
well-known differential equation 

dx^ \x-lx + ljdx'' x^ -1 " 

In this case, the Van Vleck polynomial V{x) is trivially a fc-dependent constant. Hence, there is only one 
set of zeros of p^~^'^~^^ (^Xi) with i = 1, 2, • • • , fc, satisfying the Bethe ansatz equation (16). 

The case with 7 7^ is non-trivial. According to (l)-(5), we write a differential equation corresponding 
to the Bethe ansatz equation (16) as 

d'^yix) (a P \ dyix) Vix) , , „ 



dx'^ \x — 1 x + 1 J dx x^ — 1 

with the corresponding polynomials A{x) = — 1 and B{x) = ■yx'^ + {a + B)x + a — B — j shown in 
(1). In contrast to the Heine-Sieltjes equation (1), however, in this case the polynomial B(x) is of the 
same degree as that of A{x). Therefore, the polynomials y{x) determined by (18) should be similar to 
but different from those of Heine-Sieltjes type. 

In search for polynomial solutions of (18), we write 

k 

yk{x)=J2^jX^- (19) 
Substitution of (19) into (18) yields the condition to determine the corresponding polynomial V(x) with 

V{x) = -^kx + f, (20) 

where / is an undetermined constant depending on k, together with the expansion coefficients bj, satis- 
fying the following four-term relations: 



(i(a + p+j-l) + f)bj = {j + 2){j + l)6,+2 + (j + l)(/3 - a + 7)6i+i + l{k - j + l)6,-i (21) 

with 6j = for j < — 1 or J > fc -|- 1, which is equivalent to the following eigen-equation for / with 

Fb = /b, (22) 

where 
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- a + 7) 



fc7 -(a + /3) 2(/3-a + 7) 



(A; -1)7 -2(a + /3 + l) 3(/3-Q! + 7) 12 0- 



k{k - 1) 



27 -(fc-l)(a + /3 + fc-2) fc(/3-a + 7) 



V 
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-fc(a + ^ + A;-l) j 



(23) 



and the transpose of b is related to the expansion cocfRcicnts {bj} with b-^ = (60,^1, • • • ■ bk-i, bk). It 
is clear that the matrix F should have k + 1 eigenvalues labeled as /"-^^ (C = 1, 2, • • • , fc + 1), of which 
the corresponding eigenvector b*^^^ determines the polynomial yf\x) according to (19). The number 
of solutions for / is indeed the same as that of the eigenstates of the standard two-site Bose-Hubbard 
Hamiltonian (8) when 7^0. The polynomials (19) with 7^0 determined by (21) are called the extended 
Heine-Stieltjes polynomials, which, in general, can also be obtained from the Riccati differential equation 
studied in [31-32] or from relevant Bethe ansatz equations given in [33] though these authors did not 
intend to do so. It should be noted that the polynomial approach shown above is similar to that studied 
in [34-36] for quasi-exactly solvable sextic anharmonic oscillator problem and "PT-symmetric quantum 
mechanics. Since the Bethe ansatz equations of the quasi-exactly solvable sextic anharmonic oscillator 
problem and 'PT'-symmctric quantum mechanics are different from those for the special LMG model, the 
resultants should belong to different types of extended Heine-Stieltjes polynomials. 

As is well known, the advantage of the Bethe ansatz method for Gaudin type systems lies in the 
fact that the huge matrix in the Fock subspace is greatly reduced, especially for the Gaudin-Richardson 
pairing model [22-25]. However, the non-linear Bethe ansatz equations similar to (16) are very difficult 
to be solved numerically, especially for large size systems. There are several attempts to overcome this 
difficulty. The approach via Riccati differential equation shown in [31-32] is one of them. Actually, if the 
polynomials can be derived recursively similar to (21) for the LMG model, it should be much easier to 
determine zeros from the polynomials than to solve a set of BAEs with a set of variables because there 
is only one; variable in the polynomials. In order to make this point clear, let us take a simple nontrivial 
example with k = 2 and a ^ /3 ~ •y — 1/2 corresponding to j^i = 1^2 = and t/U = 1/2. In this case, 
the three polynomials y2{x) given by (19) with 62 = 1 can easily obtained from the four-term recurrence 
relations (21) as listed in Table 1. Since we set 6fe = 1 in y{x), the coefficient bk-i must equal to negative 
sum of zeros of y{x) with bk-i = — Y^i^i Xi- Therefore, the solution corresponding to the largest bk-i is 
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that for the ground state of the system considered; that corresponding to the next largest bk-i is that of 
the first excited state; and so on. One can check that these are indeed the case as shown in Table 1. 



Table 1: The- oxt(-iulod H(nn(--Sti(-Itj(^s polvnomials ijoi-r) for tlu- LMG model with n = /3 = 7 = 1/2. 

Zeros of the polynomial xi + X2 The Polynomial 2/2(2;) / 

xi = -7.2904, X2 = -1.7379 -9.0283 12.6700 + 9.0283x + a;^ 0.5141 

= -5.7052, X2 = 0.5612 -5.1440 -3.2020 + 5.1440a; + a;^ -1.4280 

xi = -0.6036, a;2 = 0.7759 0.1723 -0.4684 - 0.1723a; + a;^ -4.0861 



According to the Stieltjes results, an electrostatic interpretation of the location of zeros of the new 
polynomial y{x) may be stated as follows. Put two positive fixed charges a/2 and /3/2 at —1 and +1 
along a real line, respectively, and allow k positive unit charges to move freely along the real line under 
such situation together with a uniform electric field with strength —7/2. Therefore, up to a constant, 
the total energy functional f7(a;i,a;2, • ■ • ,a;fe) may be written as 



C/(a;i,a;2,---,a;fc) = -|^a;, - ^^ln|a;i - 1| - ^^ln|a;i + 1| - ^ ln|a;i-a;j|. (24) 

i i i l<i<j<k 

In this case, there are A: + 1 different configurations for the position of these k charges {a;^'' ' , • • • ) x'k ^ } with 

C = 1, 2, • ■ • , + 1, corresponding to global minimums of the total energy. As proven by Stieltjes, there 

is a unique configuration with — l<a;i<---<a;/j<l when 7 = which corresponds to the zeros of the 

Jacobi polynomial P^" ^'^ When 7 > 0, however, one can verify numerically as done in [15] that 

the range of the positions of these positive unit charges tends to be along the entire half line except the two 

singular points ±1 with {xi, ■ ■ ■ ,Xk} S {—00, —1) IJ (—1, +1). Let the positions of these positive charges 

be arranged as a;i < X2 < • • • < a;fc. The above restriction requires that these positions must be within 

the intervals (—00, —1) and (—1, +1) with —00 < xi < X2 < ■ ■ ■ < Xk^ < — 1 < Xk^+i < ■ • ■ < Xk < +1. 

It follows from this that the total number of possible configurations is exactly the number of ways to put 

/fc + rn-l\ 

the k positive charges into the two intervals, which is ( ) = fc + 1 for m = 2. 

V k J 

Some explicit formulas for sums of powers of zeros of the new polynomials can easily be derived. For 
example, summing Eq. (16) over i, we have 

- /5) E ^ + (" + /?) E ^ = -^7, (25) 
while multiplying Eq. (16) by a;, and then summing over i, we get 

(" + /3)ET23T + ("-^)E:;^rT + 7E^^ = -M" + /3 + fc-i)- (26) 

2 — 1 ^ %=1 



Eqs. (25) and (26) can be combined to express 'yJ2i=i related to the eigenenergy (14) as 
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AaB 1 , a — B ^ , ^ 

-. = ^ET3^-fc(" + /3 + fc-l-7^). (27) 

Additional sum rules can be derived from the explicit form of the new polynomial yk {x) by using the 
method shown in [37-39] for other polynomials. Let the polynomial yk{x), up to a constant, be expressed 
in terms of its zeros as 

k 

ykix)=Y[{x^x^). (28) 
Then, the expansion coefficients in (19) can be explicitly written as 

k k 

bk = l, bk-i=-'^Xi, bk-2 = ^XiXj, bo = {-)'''[[xi. (29) 

1=1 i<j i=l 

Directly using the four-term relation (21) with the explicit expressions shown in (29), we have, for 
example, 

/= E ^ + ("-^-7)E^=fE^) -(fc-i)E^ + ("-/3-7)E^ (30) 

for k >2, and 

k 

f=-j^Xi-k{a + p + k-l), (31) 

i=l 

which provides a relation among the eigenvalues / in the differential equation (18) and the eigenenergies 
(14) of the special LMG model. Eq. (31) can easily be verified from the results shown in Table 1. 

Sum rules of other higher order powers and inverse powers of zeros of the new polynomials may also 
be derived from the four-term relation (21) by using explicit expressions shown in (29). Eq. (31) clearly 
shows that / = —k{a + j3 + k — 1) corresponding to the Jacobi polynomial with yk{x) = Plf~^'^~^\x) 
when 7 = 0. 



5 Summary and discussions 

In this paper, new polynomials, similar to but distinct from those of the Heine-Stieltjes type, that 
are associated with a special LMG model corresponding to the standard two-site Bose-Hubbard model 

are derived based on the Sticltjes correspondence. It follows that the eigenvalues and corresponding 
eigenstates of the special LMG model can be determined from the zeros of these polynomials. Further, 
if these polynomials can be derived recursively, it also follows that it should be much easier to determine 
zeros from the polynomials, e. g. the case considered in [40], than to solve a set of BAEs with a set 

of variables because there is only one variable in the polynomials. This conclusion applies equally well 
to other quantum many-body systems, such as atomic-molecular Bosc-Einstein condensates [41], the 
heteronuclear molecular Bose- Einstein condensate model [42], and the nuclear pairing problem, when 
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the corresponding polynomials can be obtained. It is well known that the Jacobi polynomials P^"' (x) 
are orthogonal with respect to the measure rf/x(a;) = {x — l)"{x + l)^dx on the interval x G (— 1,+1). 
Although the orthogonality of the new polynomials is not addressed, we conjuncture that, similar to the 
Lame polynomials analyzed in [43-44] , there is no sequence of the new polynomials with 7^0 orthogonal 
with respect to any measure supported on the interval (—00, 1) when 7 ^ 0. A rigorous proof of the 
latter, which requires further work, is beyond the scope of the present study. Furthermore, following 
Stieltjes early work, a one-dimensional classical electrostatic analogue corresponding to the special LMG 
model is established. This shows that any possible configuration of ciquilibrium positions of the charges in 
the electrostatic problem corresponds uniquely to one set of roots of the BAEs for the LMG model, and 
that the number of possible configurations of equilibrium positions of the charges equals exactly to the 
number of energy levels in the LMG model. Sum rules of powers and inverse powers of zeros of the new 
polynomials related to eigenenergies of the LMG model are also considered. More complicated sum rules 
and their relations may be obtained similarly. The results by extension clearly show a new link between 
Bethe ansatz type solutions of large class of quantum many-body problems and a class of polynomials 
satisiying second-order differential equations, and thus opens a new way to get solutions of these BAEs 
from zeros of these polynomials. 
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